#include "plugin.h"
#include <stdio.h>
#include <stack>
#include <QVector3D>
#include <unistd.h>

struct AlgRefineRayburst: public AlgRefine {
	static AlgRefine* creater() {
		return new AlgRefineRayburst{};
	}
	int mark;
	double radius;
	void setParams(const std::vector<ParamValue>& params) override {
		mark=params[0].getInt();
		radius=params[1].getFloat();
	}
	std::vector<ParamValue> getParams() const override {
		std::vector<ParamValue> vl;
		vl.push_back(ParamValue{mark});
		vl.push_back(ParamValue{radius});
		return vl;
	}
	AlgRefineRayburst() : AlgRefine{}, mark{0}, radius{255} { }
	void refine(Graph gr, Path* path) override {
		auto& pth=path->points;
		double res=xres<0>();
		if(res>yres<0>()) res=yres<0>();
		if(res>zres<0>()) res=zres<0>();
		res/=2;
		for(size_t i=0; i<pth.size(); i++) {
			if(pth[i]._r!=0)
				continue;
			Point p0=pth[i];
			if(i>0)
				p0=pth[i-1];
			Point p1=pth[i];
			if(i+1<pth.size())
				p1=pth[i+1];
			QVector3D dp(p0.x()-p1.x(), p0.y()-p1.y(), p0.z()-p1.z());
			dp.normalize();
			QVector3D du;
			if(fabs(dp.z())<0.57) {
				du=QVector3D(-dp.y(), dp.x(), 0).normalized();
			} else {
				du=QVector3D(-dp.z(), 0, dp.x()).normalized();
			}
			QVector3D dv=QVector3D::crossProduct(dp, du).normalized();
			std::vector<double> rvals{};
			int64_t idx_center=xyzf2index<0>(pth[i].x(), pth[i].y(), pth[i].z());
			double val_center=xfunc<0>(cube<0>(idx_center))*255;
			for(int j=0; j<16; j++) {
				double alpha=j*2*3.14159/16;
				QVector3D dw=cos(alpha)*du+sin(alpha)*dv;
				double dx=dw.x(), dy=dw.y(), dz=dw.z();
				double len=sqrt(dx*dx+dy*dy+dz*dz);
				dx/=len;
				dy/=len;
				dz/=len;
				double mdiff=1024*1024;
				double rval=0;
				int k=1;
				int val_prev=val_center;
				//double val_mean=val_center;
				for(double len=res; len<=radius; len+=res) {
					int64_t idx=xyzf2index<0>(pth[i].x()+dx*len, pth[i].y()+dy*len, pth[i].z()+dz*len);
					if(idx<0)
						break;
					double val=xfunc<0>(cube<0>(idx))*255;
					auto diff=val-val_prev;
				double val_mean=val_prev;
					if(diff*val_mean<mdiff) {
						mdiff=diff*val_mean;
						rval=len;
					}
					val_mean=(val_mean*k+val)/(k+1);
					val_prev=val;
					k++;
					if(i==0) {
						fprintf(stderr, "d%lf, ", diff);
					}
					if(val*val_mean<-mdiff)
						break;
				}
					if(i==0) {
						fprintf(stderr, "\n%lf\n", rval);
					}
				rvals.push_back(rval);
			}
			std::sort(rvals.begin(), rvals.end());
			double sum=0;
			for(int j=2; j<8; j++) {
				sum+=rvals[j];
			}
			pth[i].r(sum/6);
		}
	/*
		size_t a=0;
		size_t b=pth->path.size()-1;
		if(pth->fromEdge) {
			pth->path[a].m=pth->fromEdge.points()[pth->fromEdgePos].m;
			pth->path[a].setR(pth->fromEdge.points()[pth->fromEdgePos].r());
			a++;
		}
		if(pth->toEdge) {
			pth->path[b].m=pth->toEdge.points()[pth->toEdgePos].m;
			pth->path[b].setR(pth->toEdge.points()[pth->toEdgePos].r());
			b--;
		}
		for(size_t i=a; i<=b; i++) {
			pth->path[i].m=mark;
			pth->path[i].setR(radius);
		}
		*/
	}
};


static void astar_load(std::vector<AlgorithmInfo*>& algs, std::vector<GraphReaderInfo*>& readers) {
	auto algrefine=new AlgInfo<AlgRefine>{"Rayburst", 60,
		"Set mark and radius with Rayburst.\nIntensity values (and transfer function) of channel 1 in group 1 is used.",
		{
			ParamSpec{"Mark Bouton", "...", 0, 255, 0, ""},
			ParamSpec{"Max Radius", "...", 0.0, 1000.0, 255.0, "µm"}
		},
		AlgRefineRayburst::creater};
	algs.push_back(algrefine);
}

static void astar_unload() {
}

FNT_PLUGIN("Rayburst", astar_load, astar_unload, "Simple built-in Rayburst algorithms.");

